home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dlabrd.f < prev    next >
Text File  |  1996-07-19  |  12KB  |  292 lines
  1.       SUBROUTINE DLABRD( M, N, NB, A, LDA, D, E, TAUQ, TAUP, X, LDX, Y,
  2.      $                   LDY )
  3. *
  4. *  -- LAPACK auxiliary routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     February 29, 1992
  8. *
  9. *     .. Scalar Arguments ..
  10.       INTEGER            LDA, LDX, LDY, M, N, NB
  11. *     ..
  12. *     .. Array Arguments ..
  13.       DOUBLE PRECISION   A( LDA, * ), D( * ), E( * ), TAUP( * ),
  14.      $                   TAUQ( * ), X( LDX, * ), Y( LDY, * )
  15. *     ..
  16. *
  17. *  Purpose
  18. *  =======
  19. *
  20. *  DLABRD reduces the first NB rows and columns of a real general
  21. *  m by n matrix A to upper or lower bidiagonal form by an orthogonal
  22. *  transformation Q' * A * P, and returns the matrices X and Y which
  23. *  are needed to apply the transformation to the unreduced part of A.
  24. *
  25. *  If m >= n, A is reduced to upper bidiagonal form; if m < n, to lower
  26. *  bidiagonal form.
  27. *
  28. *  This is an auxiliary routine called by DGEBRD
  29. *
  30. *  Arguments
  31. *  =========
  32. *
  33. *  M       (input) INTEGER
  34. *          The number of rows in the matrix A.
  35. *
  36. *  N       (input) INTEGER
  37. *          The number of columns in the matrix A.
  38. *
  39. *  NB      (input) INTEGER
  40. *          The number of leading rows and columns of A to be reduced.
  41. *
  42. *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  43. *          On entry, the m by n general matrix to be reduced.
  44. *          On exit, the first NB rows and columns of the matrix are
  45. *          overwritten; the rest of the array is unchanged.
  46. *          If m >= n, elements on and below the diagonal in the first NB
  47. *            columns, with the array TAUQ, represent the orthogonal
  48. *            matrix Q as a product of elementary reflectors; and
  49. *            elements above the diagonal in the first NB rows, with the
  50. *            array TAUP, represent the orthogonal matrix P as a product
  51. *            of elementary reflectors.
  52. *          If m < n, elements below the diagonal in the first NB
  53. *            columns, with the array TAUQ, represent the orthogonal
  54. *            matrix Q as a product of elementary reflectors, and
  55. *            elements on and above the diagonal in the first NB rows,
  56. *            with the array TAUP, represent the orthogonal matrix P as
  57. *            a product of elementary reflectors.
  58. *          See Further Details.
  59. *
  60. *  LDA     (input) INTEGER
  61. *          The leading dimension of the array A.  LDA >= max(1,M).
  62. *
  63. *  D       (output) DOUBLE PRECISION array, dimension (NB)
  64. *          The diagonal elements of the first NB rows and columns of
  65. *          the reduced matrix.  D(i) = A(i,i).
  66. *
  67. *  E       (output) DOUBLE PRECISION array, dimension (NB)
  68. *          The off-diagonal elements of the first NB rows and columns of
  69. *          the reduced matrix.
  70. *
  71. *  TAUQ    (output) DOUBLE PRECISION array dimension (NB)
  72. *          The scalar factors of the elementary reflectors which
  73. *          represent the orthogonal matrix Q. See Further Details.
  74. *
  75. *  TAUP    (output) DOUBLE PRECISION array, dimension (NB)
  76. *          The scalar factors of the elementary reflectors which
  77. *          represent the orthogonal matrix P. See Further Details.
  78. *
  79. *  X       (output) DOUBLE PRECISION array, dimension (LDX,NB)
  80. *          The m-by-nb matrix X required to update the unreduced part
  81. *          of A.
  82. *
  83. *  LDX     (input) INTEGER
  84. *          The leading dimension of the array X. LDX >= M.
  85. *
  86. *  Y       (output) DOUBLE PRECISION array, dimension (LDY,NB)
  87. *          The n-by-nb matrix Y required to update the unreduced part
  88. *          of A.
  89. *
  90. *  LDY     (output) INTEGER
  91. *          The leading dimension of the array Y. LDY >= N.
  92. *
  93. *  Further Details
  94. *  ===============
  95. *
  96. *  The matrices Q and P are represented as products of elementary
  97. *  reflectors:
  98. *
  99. *     Q = H(1) H(2) . . . H(nb)  and  P = G(1) G(2) . . . G(nb)
  100. *
  101. *  Each H(i) and G(i) has the form:
  102. *
  103. *     H(i) = I - tauq * v * v'  and G(i) = I - taup * u * u'
  104. *
  105. *  where tauq and taup are real scalars, and v and u are real vectors.
  106. *
  107. *  If m >= n, v(1:i-1) = 0, v(i) = 1, and v(i:m) is stored on exit in
  108. *  A(i:m,i); u(1:i) = 0, u(i+1) = 1, and u(i+1:n) is stored on exit in
  109. *  A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).
  110. *
  111. *  If m < n, v(1:i) = 0, v(i+1) = 1, and v(i+1:m) is stored on exit in
  112. *  A(i+2:m,i); u(1:i-1) = 0, u(i) = 1, and u(i:n) is stored on exit in
  113. *  A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).
  114. *
  115. *  The elements of the vectors v and u together form the m-by-nb matrix
  116. *  V and the nb-by-n matrix U' which are needed, with X and Y, to apply
  117. *  the transformation to the unreduced part of the matrix, using a block
  118. *  update of the form:  A := A - V*Y' - X*U'.
  119. *
  120. *  The contents of A on exit are illustrated by the following examples
  121. *  with nb = 2:
  122. *
  123. *  m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n):
  124. *
  125. *    (  1   1   u1  u1  u1 )           (  1   u1  u1  u1  u1  u1 )
  126. *    (  v1  1   1   u2  u2 )           (  1   1   u2  u2  u2  u2 )
  127. *    (  v1  v2  a   a   a  )           (  v1  1   a   a   a   a  )
  128. *    (  v1  v2  a   a   a  )           (  v1  v2  a   a   a   a  )
  129. *    (  v1  v2  a   a   a  )           (  v1  v2  a   a   a   a  )
  130. *    (  v1  v2  a   a   a  )
  131. *
  132. *  where a denotes an element of the original matrix which is unchanged,
  133. *  vi denotes an element of the vector defining H(i), and ui an element
  134. *  of the vector defining G(i).
  135. *
  136. *  =====================================================================
  137. *
  138. *     .. Parameters ..
  139.       DOUBLE PRECISION   ZERO, ONE
  140.       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
  141. *     ..
  142. *     .. Local Scalars ..
  143.       INTEGER            I
  144. *     ..
  145. *     .. External Subroutines ..
  146.       EXTERNAL           DGEMV, DLARFG, DSCAL
  147. *     ..
  148. *     .. Intrinsic Functions ..
  149.       INTRINSIC          MIN
  150. *     ..
  151. *     .. Executable Statements ..
  152. *
  153. *     Quick return if possible
  154. *
  155.       IF( M.LE.0 .OR. N.LE.0 )
  156.      $   RETURN
  157. *
  158.       IF( M.GE.N ) THEN
  159. *
  160. *        Reduce to upper bidiagonal form
  161. *
  162.          DO 10 I = 1, NB
  163. *
  164. *           Update A(i:m,i)
  165. *
  166.             CALL DGEMV( 'No transpose', M-I+1, I-1, -ONE, A( I, 1 ),
  167.      $                  LDA, Y( I, 1 ), LDY, ONE, A( I, I ), 1 )
  168.             CALL DGEMV( 'No transpose', M-I+1, I-1, -ONE, X( I, 1 ),
  169.      $                  LDX, A( 1, I ), 1, ONE, A( I, I ), 1 )
  170. *
  171. *           Generate reflection Q(i) to annihilate A(i+1:m,i)
  172. *
  173.             CALL DLARFG( M-I+1, A( I, I ), A( MIN( I+1, M ), I ), 1,
  174.      $                   TAUQ( I ) )
  175.             D( I ) = A( I, I )
  176.             IF( I.LT.N ) THEN
  177.                A( I, I ) = ONE
  178. *
  179. *              Compute Y(i+1:n,i)
  180. *
  181.                CALL DGEMV( 'Transpose', M-I+1, N-I, ONE, A( I, I+1 ),
  182.      $                     LDA, A( I, I ), 1, ZERO, Y( I+1, I ), 1 )
  183.                CALL DGEMV( 'Transpose', M-I+1, I-1, ONE, A( I, 1 ), LDA,
  184.      $                     A( I, I ), 1, ZERO, Y( 1, I ), 1 )
  185.                CALL DGEMV( 'No transpose', N-I, I-1, -ONE, Y( I+1, 1 ),
  186.      $                     LDY, Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
  187.                CALL DGEMV( 'Transpose', M-I+1, I-1, ONE, X( I, 1 ), LDX,
  188.      $                     A( I, I ), 1, ZERO, Y( 1, I ), 1 )
  189.                CALL DGEMV( 'Transpose', I-1, N-I, -ONE, A( 1, I+1 ),
  190.      $                     LDA, Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
  191.                CALL DSCAL( N-I, TAUQ( I ), Y( I+1, I ), 1 )
  192. *
  193. *              Update A(i,i+1:n)
  194. *
  195.                CALL DGEMV( 'No transpose', N-I, I, -ONE, Y( I+1, 1 ),
  196.      $                     LDY, A( I, 1 ), LDA, ONE, A( I, I+1 ), LDA )
  197.                CALL DGEMV( 'Transpose', I-1, N-I, -ONE, A( 1, I+1 ),
  198.      $                     LDA, X( I, 1 ), LDX, ONE, A( I, I+1 ), LDA )
  199. *
  200. *              Generate reflection P(i) to annihilate A(i,i+2:n)
  201. *
  202.                CALL DLARFG( N-I, A( I, I+1 ), A( I, MIN( I+2, N ) ),
  203.      $                      LDA, TAUP( I ) )
  204.                E( I ) = A( I, I+1 )
  205.                A( I, I+1 ) = ONE
  206. *
  207. *              Compute X(i+1:m,i)
  208. *
  209.                CALL DGEMV( 'No transpose', M-I, N-I, ONE, A( I+1, I+1 ),
  210.      $                     LDA, A( I, I+1 ), LDA, ZERO, X( I+1, I ), 1 )
  211.                CALL DGEMV( 'Transpose', N-I, I, ONE, Y( I+1, 1 ), LDY,
  212.      $                     A( I, I+1 ), LDA, ZERO, X( 1, I ), 1 )
  213.                CALL DGEMV( 'No transpose', M-I, I, -ONE, A( I+1, 1 ),
  214.      $                     LDA, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
  215.                CALL DGEMV( 'No transpose', I-1, N-I, ONE, A( 1, I+1 ),
  216.      $                     LDA, A( I, I+1 ), LDA, ZERO, X( 1, I ), 1 )
  217.                CALL DGEMV( 'No transpose', M-I, I-1, -ONE, X( I+1, 1 ),
  218.      $                     LDX, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
  219.                CALL DSCAL( M-I, TAUP( I ), X( I+1, I ), 1 )
  220.             END IF
  221.    10    CONTINUE
  222.       ELSE
  223. *
  224. *        Reduce to lower bidiagonal form
  225. *
  226.          DO 20 I = 1, NB
  227. *
  228. *           Update A(i,i:n)
  229. *
  230.             CALL DGEMV( 'No transpose', N-I+1, I-1, -ONE, Y( I, 1 ),
  231.      $                  LDY, A( I, 1 ), LDA, ONE, A( I, I ), LDA )
  232.             CALL DGEMV( 'Transpose', I-1, N-I+1, -ONE, A( 1, I ), LDA,
  233.      $                  X( I, 1 ), LDX, ONE, A( I, I ), LDA )
  234. *
  235. *           Generate reflection P(i) to annihilate A(i,i+1:n)
  236. *
  237.             CALL DLARFG( N-I+1, A( I, I ), A( I, MIN( I+1, N ) ), LDA,
  238.      $                   TAUP( I ) )
  239.             D( I ) = A( I, I )
  240.             IF( I.LT.M ) THEN
  241.                A( I, I ) = ONE
  242. *
  243. *              Compute X(i+1:m,i)
  244. *
  245.                CALL DGEMV( 'No transpose', M-I, N-I+1, ONE, A( I+1, I ),
  246.      $                     LDA, A( I, I ), LDA, ZERO, X( I+1, I ), 1 )
  247.                CALL DGEMV( 'Transpose', N-I+1, I-1, ONE, Y( I, 1 ), LDY,
  248.      $                     A( I, I ), LDA, ZERO, X( 1, I ), 1 )
  249.                CALL DGEMV( 'No transpose', M-I, I-1, -ONE, A( I+1, 1 ),
  250.      $                     LDA, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
  251.                CALL DGEMV( 'No transpose', I-1, N-I+1, ONE, A( 1, I ),
  252.      $                     LDA, A( I, I ), LDA, ZERO, X( 1, I ), 1 )
  253.                CALL DGEMV( 'No transpose', M-I, I-1, -ONE, X( I+1, 1 ),
  254.      $                     LDX, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
  255.                CALL DSCAL( M-I, TAUP( I ), X( I+1, I ), 1 )
  256. *
  257. *              Update A(i+1:m,i)
  258. *
  259.                CALL DGEMV( 'No transpose', M-I, I-1, -ONE, A( I+1, 1 ),
  260.      $                     LDA, Y( I, 1 ), LDY, ONE, A( I+1, I ), 1 )
  261.                CALL DGEMV( 'No transpose', M-I, I, -ONE, X( I+1, 1 ),
  262.      $                     LDX, A( 1, I ), 1, ONE, A( I+1, I ), 1 )
  263. *
  264. *              Generate reflection Q(i) to annihilate A(i+2:m,i)
  265. *
  266.                CALL DLARFG( M-I, A( I+1, I ), A( MIN( I+2, M ), I ), 1,
  267.      $                      TAUQ( I ) )
  268.                E( I ) = A( I+1, I )
  269.                A( I+1, I ) = ONE
  270. *
  271. *              Compute Y(i+1:n,i)
  272. *
  273.                CALL DGEMV( 'Transpose', M-I, N-I, ONE, A( I+1, I+1 ),
  274.      $                     LDA, A( I+1, I ), 1, ZERO, Y( I+1, I ), 1 )
  275.                CALL DGEMV( 'Transpose', M-I, I-1, ONE, A( I+1, 1 ), LDA,
  276.      $                     A( I+1, I ), 1, ZERO, Y( 1, I ), 1 )
  277.                CALL DGEMV( 'No transpose', N-I, I-1, -ONE, Y( I+1, 1 ),
  278.      $                     LDY, Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
  279.                CALL DGEMV( 'Transpose', M-I, I, ONE, X( I+1, 1 ), LDX,
  280.      $                     A( I+1, I ), 1, ZERO, Y( 1, I ), 1 )
  281.                CALL DGEMV( 'Transpose', I, N-I, -ONE, A( 1, I+1 ), LDA,
  282.      $                     Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
  283.                CALL DSCAL( N-I, TAUQ( I ), Y( I+1, I ), 1 )
  284.             END IF
  285.    20    CONTINUE
  286.       END IF
  287.       RETURN
  288. *
  289. *     End of DLABRD
  290. *
  291.       END
  292.